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I. INTRODUCTION 



One of the two main avenues for current research in quantum gravity is that of quantum 
geometry (also known as loop quantum gravity). This is a theory that leads to space- 
time that is discrete at Planck-length scales, giving rise to predictions of quantized areas 
and volumes (see [1] for reviews). In recent years, these ideas have been applied to the 
study of cosmology [2, 3]. Just as minisuperspace models reduce the infinitely many degrees 
of freedom down to a finite, quantum mechanical model, the basic idea of loop quantum 
cosmology (LQC) is to use symmetry reduction to develop relatively simple models for 
various space-times. However, many of the features of LQC are similar to those in the full 
theory; in particular, the Hamiltonian constraint is constructed analogously to the full one 
[4], and as a consequence takes the form of a discrete recursion relation [5]. 

Loop quantum cosmological models are free of singularities, irrespective of the type of 
matter content [6]. This is seen as a general consequence of the quantum evolution equation, 
which is a difference equation for the wave function and does not break down where the 
classical singularity would be. At the same time, additional conditions for the wave function 
can arise as consistency conditions for solutions of the difference equation whose coefficients 
can become zero. These conditions can be interpreted as dynamical initial conditions [7]. For 
a vacuum isotropic model [9], it turns out that together with the notion of pre-classicality [7] 
a unique LQC solution is picked out. Pre-classicality requires the coefficients of the wave 
function to lack large variations over Planck-scale lengths, and thus allows the solution to 
match a semi-classical wave function far away from the singularity. Because one is dealing 
with discrete recursion relations, this means that as the parameter is increased by one, 
there should not be a great change in the value of the coefficients. (A precise definition 
of pre-classicality would require more information about observables and the physical inner 
product than is presently available [8].) 

The situation is more complicated in anisotropic models, with several independent direc- 
tions in minisuperspace along which oscillations need to be suppressed, and therefore the 
Hamiltonian constraint becomes a partial difference equation. We will look in detail at the 
Bianchi I LRS model, which has only two independent degrees of freedom and simplifies 
thanks to its separable evolution equation. As we will see, the generic sequence alternates 
between positive and negative values, thus exhibiting huge oscillations. However, with the 
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new techniques introduced here to LQC it is possible to select out special boundary values 
that give pre-classical sequences; from this set of possibilities, pre-classical wave functions 
can be built up. Because these generating function methods can be used to solve generic 
(partial) difference equations, they are useful to obtain pre-classical solutions for other LQC 
models and their Hamiltonian constraints. 

To study the relation between the quantum theory and its classical limit, it is of particular 
interest to construct solutions which in presumed semiclassical regimes describe wave pack- 
ets following classical trajectories at least approximately. This requires a sufficiently large 
set of independent pre-classical solutions such that an initial wave packet can be formed 
by superposition. That this is realized is not at all obvious. For instance, the continuum 
limit of the difference equation with the corresponding (DeWitt 1 [10]) initial condition is 
not well-posed (solutions to the Wheeler-DeWitt equation which are zero at the boundaries 
of minisuperspace do not suffice to allow reasonable initial data). Even though there are 
LQC models which are well-posed in contrast to their Wheeler-DeWitt analog [11], it is 
conceivable that in a model like that studied here the lack of well-posedness in the contin- 
uum limit precludes the existence of sufficiently many pre-classical solutions. This indeed 
turns out to be the case as will be discussed later together with possible consequences and 
interpretations. 

The loop quantized Hamiltonian constraint for the particular case of a Bianchi I LRS 
space-time used here has already been worked out [12]. Here we present more systematic 
methods of finding solutions with particular properties, which are useful not only for this 
space-time, but for any model system to be considered in loop quantum gravity. In Section 
II, we discuss in detail how to use generating functions to obtain information about a 
sequence with one parameter. This allows us to tune the magnitude and asymptotic behavior 
of the sequence by choosing the appropriate initial values. By combining two such one- 
parameter sequences, we can solve the Bianchi I LRS recursion relation, and pick appropriate 
combinations to have the desired properties at large values of the evolution parameter n. 

1 In order to have a regular quantum theory, DeWitt introduced the condition that the wave function van- 
ish at points of minisuperspace corresponding to classical singularities. Posing those initial conditions, 
however, is not enough to avoid singularities, and is replaced in loop quantum cosmology by a differ- 
ent mechanism which depends sensitively on details of a loop quantization. Still, the dynamical initial 
conditions mentioned above take a similar form and emerge automatically from loop quantum cosmology. 
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Section III generalizes the similar method to solve the full two parameter recursion relation. 
However, because of numerical limitations, the actual sequence is obtained only over a small 
extent in the parameters m and n, but still its properties can be tuned by appropriate choices 
of boundary conditions. The general solutions we obtain here are similar in form to those 
of Section II. As an application of the results we discuss the issue of pre-classical solutions 
and the semiclassical limit in the Conclusions. 

II. SEPARABLE SOLUTIONS 

The Bianchi I LRS (locally rotationally symmetric) model is the simplest anisotropic 
model, being derived from the Bianchi I model with symmetry group R 3 by imposing the 
condition of one rotational symmetry. There are thus two degrees of freedom left. In 
addition, the classical equations of motion of the model can easily be decoupled, and the 
quantum Hamiltonian separates. We will first recall the basic equations of this model from 
the appendix of [12] and then discuss its separated evolution equations. 

A. The Bianchi I LRS model 

The Bianchi I LRS model has two degrees of freedom which, in real Ashtekar variables, 
are given by two connection components (A, c) and the conjugate momenta {pa,Pc) with 
sympletic structure given by {A,pa} = and {c, p c } = 7«. Here, k = 8nG is the 
gravitational constant and 7 the Barbero-Immirzi parameter of loop quantum gravity. The 
momenta (pa, p c ) are components of an invariant densitized triad which determine the scale 
factors a j of a Bianchi I metric by a x = -\/bc| = 02, 03 = Pa/ \/\Pc\- Thus, the metric is 
ds 2 = \p c \(dx 2 + dy 2 ) + p 2 A /\Pc\dz 2 in Cartesian coordinates and is degenerate for p A = or 



The behavior of (A, c) and (pa,Pc) as functions of time is determined by the Hamiltonian 
constraint 



which is proportional to — ApA~ 2cp c , leading to decoupled Hamiltonian equations of motion 
solved by pa oc y/r, p c oc r 2 , A oc l/y/r, c oc l/r 2 . After introducing a new time coordinate 
t := r 3 / 2 , we obtain pa oc t 1 / 3 and p c oc t 4 / 3 and thus scale factors a\ = a 2 oc t 2 / 3 , a 3 oc £ -1 / 3 



Pc = 0. 



H = 



A(2cp c + Ap A ) 



(1) 
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whose exponents ot\ = a 2 = 2/3 and 0:3 = —1/3 indeed solve the Kasner conditions for 
Bianchi I solutions, ai = 1 = J2 T aj in the special LRS case ol\ = a 2 - 

After quantizing the model with loop techniques [12], the triad components pa and p c 
become basic operators with discrete spectra, pA\fn,n) = j , y£ 2 pm\'m,n) and p c \m,n) = 
^ , ~f£'^n\m,n). The wave function is thus supported on a discrete minsuperspace, and as a 
solution to the quantized Hamiltonian constraint subject to a difference equation, 

,n—l tm— l,n+l ^m— l,n— 1 ) = 0, (2) 

where t m ^ n are the coefficients of the wave function (rescaled by a factor of the world volume 
of each basis state). Also we have 

f n = 

d{n) = < i i (3) 

and 

, N f m = 

4H= , (4) 
U m > 1 

m — 

For the remainder of the paper, the parameter n will act as a "time" parameter. Defining 
t m ,n = t m +i, n — t m -i,n (jn > 1), the recursion relation simplifies to 

d(n)[t m+hn - t m - hn ] + 2d 2 (m)[t m: n +1 - t m , n _i] = 0. (5) 

Note that for typographical simplicity, here and throughout the rest of the paper the notation 
t mtn and t mtn has been reversed from their use in [12]. Since t m , n includes the volume of each 
basis state, it follows that i m>n must vanish at the boundaries, to,n = = t m p. Because of 
its definition in terms of t m „, the values t ,n are freely specifiable; the boundary condition 
to,n = is then used in computing t from t via i m +i, n = t m ,n + tm-i,n- We still must have 
t m ,o = 0, which will act as boundary conditions for our sequence. 

The continuum limit of the difference equation is obtained at large triad components (or 
in the 7 — > limit [13]), which implies m, n ^> 1 such that the difference operators can be 
approximated by differentials. One thus arrives at the Wheeler-DeWitt equation 

W^(PA, Pc ) + 2 P 2 1 ^4(p A , Pe ) = , (6) 
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where is the continuous function that interpolates the discrete t m , n for large m, n. Note 
that the difference equation (2) for t m „ is of order four, higher than that of the Wheeler- 
DeWitt equation, which is always second order. This is a consequence of the loop quan- 
tization and cannot be avoided. It implies that there are more independent solutions to 
the discrete equation than the Wheeler-DeWitt equation has; in fact, not all the discrete 
solutions have a continuum limit. Those which do have such a limit are the pre-classical so- 
lutions which do not change rapidly if the discrete labels m and n are increased. It is known 
that, when boundary conditions are not taken into account, there are always sufficiently 
many pre-classical solutions for the semiclassical limit to be achieved [14]. In practice, find- 
ing those solutions is usually difficult when also boundary conditions have to be taken into 
account. 

Solutions of the Wheeler-DeWitt equation can easily be studied after introducing 
i>(pA,Pc) '■= dip/dpA and separating the resulting differential equation PAdip/dpA + 
4p c dip/dp c = 0. Writing ip(pA,Pc) = &{pa)P(j>c) we obtain the conditions pa®.' = \a and 
p c /3' = — 4A/3 solved by a(pA) oc p\, (3(p c ) oc p~ AX - Thus, any non-zero solution ip diverges 
either at the boundary pa = or at the classical singularity p c = unless A = 0. For 
the original wave function ip = J ipdpA this implies that solutions regular at the boundary 
are available only for — 1 < A < 0. DeWitt's condition i]j(pa,0) = ^>(0,p c ) = in this 
model, which says that there is zero probability of finding the universe at a singularity, is 
thus ill-posed since not enough separable solutions for constructing given initial values are 
available [2]. 

The difference equation (5) for t mjn can be solved similarly by assuming that the sequence 
t m , n is separable into two one-parameter sequences a m and b n , i.e. t m , n = a m b n . This is 
possible if a m and b n satisfy 

2A 

m 

b n +i ~ K-i = -\d(n)b n , (7b) 

where A is the separation parameter now for the sequences. Any choice of these two sequences 
will solve the original recursion relation; because the a m and b n relations are second-order, 
the first two values ao and a±, for example, are enough to describe the rest of the sequence 
a m for a particular A. However, as noted before the order is higher than that of the 1st 
order separated Wheeler-DeWitt equations and there are additional solutions with large 
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FIG. 1: Oscillatory behavior of a generic solution to the recursion relation (7a) for the sequence 
flm(A = —1). The initial values are ao = a± = 1. 

oscillations. For negative A, for instance, picking any two random values for the boundary 
data will generically give sequences that alternate between positive and negative values with 
magnitudes that increase as the parameter (m or n) increases. An example of such a generic 
choice is shown in Figure 1. Because of this alternation, the wave function will not be 
smooth at large n and therefore not be pre-classical. 

That solutions generically have alternating behavior can be understood as follows: For 
large m, the right hand side of (7a) is usually small compared to the left hand side such 
that a m+2 ~ a m . The relation between neighboring values, a m+ i and a m , is determined by 
the initial values a$ and a\. Together they determine a-i = 2Xa\ + a®. For negative A, and 
a = 0, say, 0,2 has the opposite sign from a\ which translates to a m having the opposite 
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sign from a m+ \ for large m and thus alternating behavior. If ao 7^ 0, a 2 may have the same 
sign as a±, but still generically oscillations will set in later. Only for positive A is it easy to 
suppress the oscillations. However, since A enters the equation for b n with the sign reversed, 
now the 6-sequence will develop alternating behavior such that the full solution alternates 
generically. 

We will develop techniques, which can be used analytically or numerically, to check 
whether it is possible to choose special initial values such that oscillations are suppressed. 
This will fix the relation between initial values and thus reduce the amount of freedom in 
pre-classical solutions. In particular, values like ao may not be allowed to vanish. While this 
does not pose a problem for the a-sequence, for the 6-sequence we have a different initial 
value problem, requiring b = as a consequence of t mj0 = 0. Thus, for the b n the only way 
to suppress oscillations is by restricting to A < 0. 



B. Generating functions 

Our strategy for choosing the initial values of the sequence a m will be to work with a 
generating function 2 . This is a function of one variable x such that the value a m is the 
coefficient of x m in a Taylor expansion, i.e. 

00 

F(x) = Y,a m x m . (8) 

m=0 

This will allow us to go from a recursion relation for a m to a differential equation for F(x). 
To see how this works, we look at how derivatives act on the function x k F(x), for a fixed k. 
We have 

= X> + k)a m x m+k -\ (9) 

so we can see that linear functions of m in the recurrence relation become linear derivatives 
of the generating function. This makes it easier to work with a m instead of b n because ^(m), 
unlike d(n), is polynomial and the map between recursion relation and differential equation 
is more obvious. Techniques for studying the values of the b n sequence will be discussed 
later. 

2 For a good reference on generating functions, see Wilf [15]; the related notion of Z-transforms is covered 
in, e.g., Oppcnhcim and Schafcr [16]. 
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We start with the recursion relation (7a) for a m , multiply it by mx m 1 and sum over all 
values of m for which the relation is valid (i.e. m > 0); then the sum is written as 3 



y~][( m + l ) a m+2 - 2Aa m+ i - (m + l)a m ]x m = 0. 



(10) 



m=0 



By mapping instances of m into derivatives, we arrive at the differential equation for F(x) 
given by 



or 



d 

dx 
d 

dx 



F(x) - a 



x 



— xF(x) 



-2X F{x) - a °=0 



x 



l-x< 



X 



F(x) 



2A £M + ao i±a£ = o. 



X 



r 



(11) 



(12) 



As written, this equation has singularities at x — — 1, and 1, which makes it problematic to 
expand around x = in the Taylor series expansion (8) of the generating function; however, 
we are interested not in the solution itself, but in the relation between the coefficients in its 
series expansion. Because of this, there is some freedom to find an equation more tractable 
to analysis, so it is natural to define a new function 

F(x) - a 



G(x) = 



x 



(13) 



Substituting this in to the relation for F(x) gives a simpler differential equation for G(x), 



d 

dx 



(l-x 2 )G(x) 



2XG(x) = oq. 



(14) 



By going from F(x) to G(x), we have reduced the number of singularities by one, but the 
question now is to relate this new generating function to the sequence a m . If we take 



G(x) = a mX m , 



(15) 



m=0 



then the mapping (13) between the generating functions F(x) and G(x) implies that the 
two sequences a m and a m are related by 



a m = a m -i, for m > 1. 



(16) 



3 We start the sum at m = 0, so that all coefficients are of the form a m +k, instead of beginning at m = 1 
(with a a m _i term) to avoid placing extra relations on the coefficients a m . This has the effect of shifting 
m — > m + 1. 
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Thus, simplifying the differential equation has resulted in a "shift" upward in the sequences. 
Once we know the two initial values of the a m sequence, ao(= °i) an d «i(= a 2 ), then we 
can find those of the a m sequence by using 

«i — a Q = 2Xa . (17) 

which comes from the recursion relation (7a) for a m in the particular case m — 1. Thus, 
the relation between a and a± implies one between ao and a±. The information in this 
differential equation is equivalent to that in the recurrence relation; any series solution of 
the differential equation (14) will solve the relation (7a) for a m , after shifting a m back to 
a m via (16). 

C. Asymptotic behavior 

The question now is how to avoid alternating oscillatory behavior in the sequence. We 
look at the function (1 — x)G(x), which is a function that generates the differences between 
adjacent values of a m , that is, 

oo 

(1 - x)G(x) = a + y^(a: m+ i - a m )x m+1 . 

m=0 

If this function has no singularities, then the sequence a m will converge to a finite value 
without oscillation. Because of the (1 — x 2 ) factor appearing in the differential equation 
(14) for G(x), the function will in general have poles at x — ±1. Thus the function is 
regular for \x\ < 1, and if there is no singularity at x — 1, we have (1 — x)G(x)\ x= i = 
a o+Em=o( a ™+ 1_ a ™) = nm m-+oo a m such that the sequence a m has a finite limit. Moreover, 
if there is not a pole at x — —1, lim„ woo (o: m+1 — a m ) = 0, so alternating oscillations are 
suppressed. This avoids situations such as that seen, for example, in the Taylor expansion 
of (1 + x)' 1 , where the coefficients of x m alternate sign: 

_L_ = i _ x + x 2_ x 3 + ... _ / lg x 
1 + x 

Obviously, (1 — x)G(x) will have the same behavior at x — — 1 that G(x) does; we avoid 
oscillations by requiring no singularity in the function G(x) at x — —1. From this, we get a 
relation between the two initial values a and a±. The story is slightly different at x — 1; 



10 



G(x) may blow up there, but (1 — x)G(x) can be finite if, e.g., G(x) has a simple pole at 
x — 1. Two examples of this are seen in the expansions for (1 — x) -1 , 

1 



1 -x 

and ln(l — x), 



1 + X + X A + £ d + 



— ln(l — x) = x + -x 2 + -a; 3 + • • • . 

In the first, the coefficients of the sequence maintain a constant value; in the second, they 
converge to zero. We will see that the latter is exactly the behavior found generically in the 
solutions G(x) for negative A. 

With this in mind, we take as a particular example the choice A = — 1; solving our 
differential equation for the generating function G(x) given in (14) with the condition G(0) = 
a , we obtain the solution 

_ a - (2a + ai)x - (4a + 2c*i) ln(l - x) 
(1 + x) z 

In obtaining this solution from the differential equation (14) for G(x), we have used the 
relation (17) between a and the initial values a and a±. For generic values a and a±, this 
function will have singularities at x — ±1; to ensure that the singularity at x — — 1 does not 
give rise to oscillatory behavior at large m, we require that 

lim [(1 -x)(l + x) 2 G(x)\ = 0. (20) 

x^ — l 

When we solve this relation, we find that 

/41n2-3\ 

which already implies that (1 — x)G(x) is regular at x — — 1. Notice that the x = 1 
singularity remains because of the ln(l — x) term; however, (1 — x) ln(l — x) is zero at 
x = 1, so (1 — x)G(x) is regular. As discussed above, this is a "good" singularity where 
the coefficients of the Taylor series go as 1/m; this behavior in a m is seen in Figure 2. We 
remark here, however, that only for negative A do we get the type of behavior seen above; 
when A > 0, the pole at x — 1 is of higher order. 

For arbitrary A we first solve the homogeneous equation (14) for a = by G(x) = 
c(l + x) A_1 (l — x) _A_1 . By varying the constant c we obtain the general solution to the 
inhomogeneous equation as 

G(x) = c(x){l+x) x -\l-x)- x - 1 (21) 
11 
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FIG. 2: Values of the sequence a m (A = —1) where the initial values are chosen to satisfy (4 In 2 — 
3)ao = (2 In 2 — \)ol\. This ensures the sequence converges to a constant value at large m. 

with 

c(x) = a J (ytO dt = C °~ X^T (1 + _ A ' ~ A; 2 ~ A; (1 + x )l 2 ) ( 22 ) 

for A 7^ 1 in terms of the hypergeometric function 2^1- (If A = 1, the equation can be 
integrated in a manner similar to the A = — 1 case.) This gives 

G(x) = c (l + x) A - x (l - x)-^ 1 - -^-(1 - x)- x -\F l {\ - A, -A; 2 - A; (1 + x)/2) (23) 

A — 1 

where only the first term determines the singularity structure at x = — 1 since the hyperge- 
ometric function 2^1(0, b; c; z) is regular at z = taking the value ^(a, b; c; 0) = 1 for all 
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a, b, c. Thus, the singularity at x — — 1 can always be removed by choosing Co = 0. Since 

ai = a = G(0) = c - 2 A a /(A - 1)^(1 - A, -A; 2 - A; 1/2) 
= c + a - Aa (^(l/2 - A/2) - ^(1 - A/2)) 

with the digamma function if)(z) = dT(z)/dz, this translates to the condition 

ai = oq(1 - A(^(l/2 - A/2) - ^(1 - A/2)) . (24) 

This expression is finite for all A which are not positive integers since the digamma function 
is analytic except for simple poles at —z G N. For A = — 1, for instance, we can use 
ip(l) — ip(S/2) = 2 In 2 — 2 and re-obtain the special case studied before. 

At x — 1, 2^1 ( a ) b; c; (1 + x)/2) has a branch point which is logarithmic for c — a — b e Z 
or c — a — 6 ^ Q. Thus, (1 — x)G(x) always has a singularity at x — 1, which for positive A 
is enhanced by the factor (1 — x)~ x_1 . Thus, for A > the sequence a m is unbounded. 

Now we turn to putting the a m and b n sequences back together. For each A, a pre-classical 
sequence t mj „ - when the boundary values have been appropriately tuned - is completely 
characterized by the choice of A, with an additional scaling factor. Because b = by 
the boundary conditions on the sequence, changing b\ will simply scale the magnitude of 
the b n values. Similarly, pre-classicality requires a relation between ao and a±, so one is 
fixed by the choice of the other, while varying this choice will again only scale the a m 
sequence. Since increasing bi by a constant factor and decreasing both ao and a\ by the 
same factor (to preserve their relationship) leaves the £ m>n sequence constant, only one of 
these is independent. Thus, a general pre-classical solution t m>n can be written as a sum 
^ A<0 c(A)t mi „(A), with the coefficient c(A) used to scale the individual terms of the sum. 
We can use linear combinations of these to get a generic solution. As an example of this, we 
require that the sequence i m>n match a Gaussian wave packet in the parameter m, at large 
values of n (see Figures 3 and 4). One difficulty in doing this, however, is apparent when 
looking at the a m graph, Figure 2; because the sequence becomes essentially constant for 
large m, it is difficult to get the necessary discrimination in m for the wave packet 4 . 

4 Recall that the sequence t TOi „ was the "m derivative" of the original sequence f m , n , so one might wonder 
if these conclusions change when t m ^ n is used. In fact, they remain much the same - since a m gives the 
change in t m ^ n with m, its decline to zero means that t m . n is basically constant for large m. Thus, again 
only at small m is the spatial discrimination possible to build a wave function with desired properties at 
large n. 
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FIG. 3: Values of the sequence t m ^ n where a linear combination of solutions (imbn were chosen such 
that the wave function evolves into a Gaussian wave packet at large n. 

III. GENERAL SOLUTIONS 

In the previous section, we looked at separating out the m and n dependence into two 
distinct sequences, then using generating function methods to analyze the a m sequence. 
Because this has reduced the problem to finding the right initial values and a 1; there is 
much to be said for the simplicity of this method. However, we can ask whether the full 
recursion relation (5) is amenable to study by generating functions. We shall see that this 
is feasible, and provides a method to study other spacetimes where the recursion relation 
may not be separable (which includes the Schwarzschild black hole interior [17]). 

Starting again with the general difference equation (5), we can work with a generating 
function F(x, y), which describes the sequence. Because occurrences of n translate to deriva- 
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FIG. 4: Close-up of the sequence given in Figure 3, showing the Gaussian wave packet at large 
values of n. 

tives d/dy of the generating function, the function d(n) on the surface presents a problem - 
there would be a non-polynomial function of derivatives acting on the generating function. 
To deal with this, we write the sequence £ OT)Tl as 



m,n 

k=0 



y^m.n' (25) 



where the tm,n satisfy 

J_r f (o) (o) , 2_r(0) _ (0) ,_ n (26 ) 

2 rm+l,li b m-l,n\ ' rm,n+l °m,n—l\ ~~ v \4\}<t>) 

J_r £ (fe) _ Ah) , 2^r.(*) _ ,(k) i _ / J_ _ A / ,(k-l) _ Ak-l) s f26b) 

2„L ( 'm+l,n ( 'm-l,nJ T L t m,n+1 ''m.n-lJ ~~ I 2n / m+1 ' n m-l,n) ■ y^uuj 

Here, we are using the fact that d(n) ~ l/2n to write the sequence t mjTl as a perturbation 
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series. In other words, we choose the leading order sequence tm,n to satisfy a recursion rela- 
tion with coefficients polynomial in the parameters m and n, such that a partial differential 
equation can be found that is equivalent to the t^m,n relation (26a). Boundary or initial 
values would then be set such that is pre-classical, i.e. by removing singularities of the 
generating function. For higher terms t^ k \ k > 1, we can choose vanishing boundary and 
initial values. Each such contribution by itself would not be pre-classical, but its oscillations 
would contribute only small perturbations to . Indeed, because of the small difference 
between l/2n and d(n), the right hand side of the equations for tm]n, k > is suppressed by 
a factor of n~ 3 compared to the left hand side. It is thus consistent with pre-classicality to 
use the relation for in order to fix boundary values as before, and then solve the relations 
for higher with boundary and initial values zero. With a vanishing right hand side the 
solution would be tm] n = for k > 1, and the small right hand side will lead to small values 
for the higher terms which become negligible compared to the tm] n values. Because of this, 
for the rest of the paper we focus on solving for the sequence tm,n', obtaining finer accuracy 
for small m, n would be possible using higher tm] n 

Working with the tm,n sequence, we look for a partial differential equation for the gener- 
ating function 

oo oo 

%^EEftA n (27) 

m=0 n=0 

Multiplying the recursion relation (26a) for tm,n by (l/2)mnx m ~ 1 y n ~ 1 and summing over all 
values it is valid gives 



X TV . 



EE 



x m y n = 0, (28) 



m=0 n=0 

and we arrive at an equation for F(x,y) which is similar in form to the separable case, i.e. 
there are singularities at x, y — 0. To avoid this, we define a new function G(x,y) as 



G(x,y) = — 
xy 



F(x,y) -J2Cv H 



[F(x,y)-yA(y]] 



xy 



My) = J2Cy n ' 1 (3°) 



,(o) 

n=l 

where 

oo 
n=l 

is related to the generating function of the m = boundary for the original tm,n sequence 
(the definition of A(y) chosen here is to simplify the form of our eventual differential equation 
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for G(x,y)). This sum begins at n — 1, since the boundary conditions imply r ' = 0. This 
gives a "shifted" sequence, since 

F(x,y) = xyG(x,y) + yA(y). (31) 

Note the similarity between this relation between F(x,y) and G(x,y), and the analogous 
one between F(x) and G{x) in the separable case. Again, the coefficients are shifted by one, 
now in each of the parameters; if 

oo oo 

G(x,y) = J2J2P^ m y n (32) 

m=0 n=0 

then 

fffn=Pm-l,n-l. (33) 

This gives us the PDE 

" x 2 )G(x,y)] + |-[(1 - y 2 )G(a;,y)] = A(y). (34) 

Again, we will have singularities at both x = ±1 and y — ±1. 

Our discussion of the singularities in Section II for separable solutions carries over here 
similarly for functions of two variables. First, we want the functions (1 — x)G(x,y) and 
(1 — y)G(x, y) both to be finite for all values of x and y, respectively, to ensure that at large 
m or n, the sequence is either constant or smoothly changing. Putting these together means 
that (1— x)(l— y)G(x, y) has to be finite for all x, y. As with the separable solutions, the only 
problem with this might be at x = — 1 and y — — 1. One condition for (1 — x)(l — y)G(x, y) 
to be finite along these two lines in the x — y plane is that the function 

H(x,y) = (l-x 2 )(l-y 2 )G(x,y) (35) 

satisfies H(—l,y) = H(x, —1) = 0. When we solve for G(x,y) in terms of H(x ) y) in the 
above relation (35), we find the new equation 

i (1 _ x ^_EiMl + (1 _ y ^_mMl = (1 _ X , ){1 _ My) , (36) 

So, given our boundary conditions on H(x, y) mentioned above, and any choice for A(y), we 
can solve this equation numerically. 

This is exactly how we would go about finding the sequence if there was a desired sequence 
of values t ,n, which is what A(y) specifies. However, more likely one wants a particular 
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limiting case for the sequence t m:Jl for large n, for example. Here, the properties of the 
functions G(x, y) and H(x, y) come into play. Because it is assumed that (l—x)(l—y)G(x, y) 
is finite, then G(x, y) itself can have simple poles at x — 1 and y — 1. For the cases mentioned 
in Section II - the functions (1 — x)^ 1 and ln(x — 1) - the coefficients in the Taylor series 
converge to constant values (unity in the first sequence, zero in the second). So, if one looks 
at the expansion in y, the coefficients of G(x, y) will assume some constant profile in x for 
large powers y n ; this profile will be determined by, e.g., the residue of the function G(x,y) 
at a simple pole y — 1. Up to a constant factor, this is simply the function H(x, 1); along 
the lines x — 1 and y — 1, the function H(x,y) gives us information about the asymptotic 
behavior of the sequence t m ri . In fact, H(x,l) is simply the generating function reflecting 
how t m ^ n behaves at large n, and similarly for H(l,y). Thus, to assemble the desired profile 
at, say, large n, one requires that H(x, 1) be a certain function, treating the y n coefficients 
of the source term A(y) as unknowns to be solved for. This can be done either numerically 
or by using pseudo-spectral methods; the form of the partial differential equation makes it 
particularly amenable to the use of Chebyshev polynomials. An example of this is shown 
in Figure 5; notice that this has the same basic shape as the sequence assembled out of the 
separable solutions in Section II. 



IV. CONCLUSIONS 



Loop quantum cosmology has discrete evolution equations of higher order than the 
Wheeler-DeWitt differential equation which is obtained in a continuum limit. There are 
thus additional independent solutions which display an alternating behavior. Those solu- 
tions are not necessarily unphysical and can in fact contain important information about 
quantum effects. Nevertheless, for studying the semiclassical limit of those models it is 
important to have control on the subset of solutions for which the small-scale oscillations 
are suppressed. In particular, one needs to find selection criteria for this subset, usually by 
specifying initial or boundary values. Choosing those special values, if possible at all, is thus 
not done to achieve a particular physical consequence and does not amount to fine-tuning. 
As demonstrated here for the Bianchi I LRS model, appropriate solutions can be extracted 
analytically by employing generating function techniques. It is emphasized here that these 
methods will work with other partial difference equations, so they can be applied to any 
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FIG. 5: Values of the sequence T m n obtained from the Taylor expansion of the generating function 
G(x,y). Notice the similarity in shape with Figure 3 obtained by combining the one-parameter 
sequences a m and b n . 

LQC model. 

For this particular model we can apply the results to discuss the initial value problem. 
In general, one would need separated solutions for all values of the separation parameter 
in order to construct a given initial wave packet at large n, far away from the classical 
singularity at n = 0. The Wheeler-DeWitt equation with DeWitt's initial conditions does 
not allow sufficiently many solutions and thus presents an ill-posed initial value problem. 
For the difference equation of loop quantum cosmology, the initial value problem is well- 
posed in the sense that there are always non-trivial solutions. However, in the semiclassical 
limit one also has to find enough solutions which do not oscillate strongly at small scales. 
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One may then speculate that analogously to the ill-posed Wheeler-DeWitt problem there 
are not sufficiently many pre-classical solutions. Another system which is ill-posed from the 
Wheeler-DeWitt point of view, an isotropic model with a free, massless scalar has been 
studied in [11]. There it turned out that the loop problem is well-posed and still allows 
enough pre-classical solutions. However, the situation there was different in that the matter 
term was responsible for the ill-posedness of the continuum formulation. The Bianchi I LRS 
model thus presents a qualitatively different system to study well-posedness, which can be 
done in detail with the methods developed in this paper. 

What we have seen is that for negative separation parameter A we can find pre-classical 
solutions. The problems of the Wheeler-DeWitt quantization at the boundary = do not 
occur in the discrete formulation and pre-classicality puts strong restrictions on the sequence 
a m fixing its boundary values. The sequence b n has only small oscillations for A < such 
that one can choose a pre-classical solution. 

For positive A, on the other hand, the situation is different. Now, the a m sequence is not 
problematic, and for b n one would have to choose special initial values in order to guarantee 
pre-classicality. However, since bo = is required by the difference equation, this option 
is not available and b n will always be alternating for positive A. Thus, b n does not have 
a continuum limit but (— l) n b n does. However, the corresponding continuum limit of the 
difference equation acting on (— l) n a m b n has an additional relative minus sign from the (— l) n 
such that one obtains a different Wheeler-DeWitt equation: Introducing T m>n := (— l) n t mj „ 
into the difference equation (5) for t mjTl , we obtain 



where the sign of the second term is switched. This is an example of different continuum 
limits obtained from one and the same difference equation, which, from the point of view 
of the discrete formulation, can be interpreted as a duality between different continuum 
formulations [5]. 

In order to construct semiclassical wave packets at large n, one can try to use only negative 
A in order to avoid oscillations. However, the resulting solutions are mostly concentrated at 
small values of m since the pre-classical separable solutions a m decrease with m. One can 
also see that using only negative A is not enough by observing that the Wheeler-DeWitt 



(-l) n d(n)(T m+1 , n 
(-l) n [d(n)(T m+1> 
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solutions p\ for positive integer A are just what one needs for a Taylor expansion from 
which one could construct any analytical initial wave packet. Using these A in the discrete 
case requires to have alternating b n , which implies that such wave packets do not follow the 
classical trajectory since they are solutions to the Wheeler-DeWitt equation with a sign flip. 

Thus, in this model the discrete analog of De Witt's initial condition may be problematic. 
It is possible to avoid DeWitt's initial condition by using a symmetric version of the con- 
straint (which is non-singular if the constraint is symmetrized after multiplying with sgn(n) 
[9]). In this way detailed studies of models can teach lessons for the full theory since the 
same ordering should be used in all cases. With this symmetric ordering, 6 would be free 
such that it can be fixed as done for the a m sequence if one requires a pre-classical solution. 
However, due to the symmetrization of the constraint operator the difference equation will 
change and be more complicated; moreover, one loses conditions on the wave function. Even 
with the non-symmetric constraint used here, it follows from the results of this paper that, 
with the condition of pre-classicality in addition to dynamical initial conditions, the wave 
function will not be unique. 

Alternatively, one can take the lack of pre-classical solutions seriously for this model. If 
we only look at large values of m and n, i.e. only at the semiclassical regime, there is no 
problem at all and we can construct any initial state we like. Most of those states would not 
solve the initial conditions imposed by the quantum constraint at the classical singularity 
and thus have to be discarded. At this point, however, we already use information from the 
quantum theory: the behavior of the wave function at the classical singularity is invoked. 
Thus, the properties of wave functions discussed here can be interpreted as a quantum effect. 
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